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In this article, we have proposed an epidemic model by using probability cellular automata theory. 
The essential mathematical features are analyzed with the help of stability theory. We have given an 
alternative modelling approach for the spatiotemporal system which is more realistic and satisfactory 
from the practical point of view. A discrete and spatiotemporal approach are shown by using cellular 
automata theory. It is interesting to note that both size of the endemic equilibrium and density 
of the individual increase with the increasing of the neighborhood size and infection rate, but the 
infections decrease with the increasing of the recovery rate. The stability of the system around 
the positive interior equilibrium have been shown by using suitable Lyapunov function. Finally 
experimental data simulation for SARS disease in China and a brief discussion conclude the paper. 

PACS numbers: 05.50.+q, 87.23.Cc, 87.18.Hf, 89.75.Fb 
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I. INTRODUCTION 



The foundation of modern mathematical epidemiology 
based on the compartment models were laid in the early 
20th century. From the middle 20th century, mathemat- 
ical epidemiology has been growing exponentially. Vari- 
ous mathematical models have been proposed and ana- 
lyzed. Although epidemic models are studied too much 
from that time but little attention has been paid on the 
localized processes of pathogen transmission between sus- 
ceptible, exposed and infected individuals. In the epi- 
demic model, generally, it is assumed that the population 
is a continuous entity and as a result, it is often neglected 
that the population are composed of single interacting 
individuals. In fact, the spread of disease must be gov- 
erned by the localized process [1] . The cellular automata 
(CA) or lattice gas cellular automata (LGCA) is a dis- 
crete approach to the time, spatial and state. It shows 
the features of the epidemic model clearly with mathe- 
matical analysis and numerical simulations in a very sim- 
ple way [2]. In a class of spatial epidemic models, it is 
supposed that the lattice of habitable sites, where no 
more than one individual can occupy any particular site. 
Spatially structured pathogen transmission says that the 
probability of susceptible host acquires some disease de- 
pends more on the local density of the exposed and in- 
fected host than that of the global [3]. However, in most 
of these models, local host density is treated as an equiv- 
alent to the global host density, since every lattice site 
is occupied with equal probability. The following spatial 
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heterogeneity is generated by exogenously, H. The va- 
rieties of the local host densities are very important to 
understand the role of spatial heterogeneity in the plants 
and animals communities [4, 5]. 

Harcourt [6] and Kobayashi [7] have studied an epi- 
demic model on the Citrus Variegated Chlorosis (CVC) 
by considering spatial dispersion and obtained a interest- 
ing negative binomial distribution. Duryca [8] has stud- 
ied an SIS epidemic model with varying pupolation den- 
sity and obtained the condition for which the pathogen 
extinct. Similar results have been shown by MoUision [3] 
and Ahmed and Agiza [9]. Keeling [10] have studied the 
effects of local spatial structure on epidemiological inva- 
sions in an epidemic model. Here we have considered an 
SEIS epidemic model with infectious force in the exposed 
and infected period, like as in [11, 12]. In this SEIS model 
we have considered the spatial dispersion and obtained 
the MF transmission threshold for epidemic persistence, 
TZc, by using probability cellular automata. We have also 
shown the nature of solution of the system around the 
trivial equilibrium (0,0) by using Lyapunov function. Fi- 
nally numerical simulations have been performed by the 
help of experimental data for the SARS disease in China. 



II. ANALYSIS FOR THE MODEL 

In order to get the same neighborhood for each site of 
the boundary, we treat the plane as a two-dimensional 
torus with J sites, where J ^ 1. The constant J is the 
carrying capacity of the environment, which is usually 
determined by the availability of the resources. In more 
general cases, the neighborhood for each site k are shown 
in Figure 1 and Figure 2. The Figure 3 shows the same 
fact with a special black point neighborhood for site k. 
Let us suppose that the individual dispersion in the sites 
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of the two-dimensional torus for the global state set is 
dented by Z{t) at time t. Then we have 

Zit) = {zi{t),Z2{t),--- ,Zj(t)}. 

In the above global set. the Zk{t) denotes the k-th set 
of local states at time t. Each site of J is empty 
or occupied by at most one individual. So, in the 
epidemic each site state belongs to the set S, where 
S = {susceptible, exposed, infected, empty}. Hence, each 
Zk{t) belongs to a finite set of elementary states. The 
recoveree from pathogen infection docs not confer immu- 
nity, so the recovered individuals will again be treated as 
susceptible individual. 

Let us assume that there exist N individual habitable 
sites ( or cells) at time t and is denoted by N{t). Each 
of the N site contains a single host, where no more than 
one individual can occupy any particular site. Then there 
are J — N{t) number of empty sites at time t, where the 
empty site represents no individual host is there. We 
denote H{t) as the global density for the individuals at 
time t. Since J is the maximum number of individual for 
the system. So we have 
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(1) 



Note that in the Eq. (1), J > N{t), therefore wc get 
H{t) < 1. In this model, we assume that the natural 
birth rate equal to it's natural death rate and d is the 
mortality rate due to the infection. 



FIG. 1: Von Neumann Neighborhood. 



FIG. 2: Moore Neighborhood. 



FIG. 3: r = 1/2. 
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FIG. 4: Schematic flow diagram 

of transmission for SEIS model. 
Here the 5 represents the maximum number individual 
of in neighborhood for site k, then the number of sites in 
5k \s 5 — \5k\ — 1, where the 5k is the k-th neighborhood. 
There is no variation in size among the sites. The local 
transition rule governs the infection transmission about 
any site k. Let the mode of transmission of the infection 
within 5k like as follows: (i) if the k-th site is non-empty, 
then we can assume that a exposed or infected k-th site 
transmits the pathogen to a susceptible individual within 
5k and (ii) if the k-th site is not infected or exposed, 
then it acquires the pathogen ( infection) from a exposed 
or infected individual within 5k- In our model, let us 
assume that 5k is symmetric. We also assume that the 
site k (or k-th site) with it's center is like as it is shown 
in Figure 1 and Figure 2. There are the same transmit 
rule for all 5k- When the 5k is asymmetric, a host at site 
k is infected from a sets differing from the set to which 
the pathogen can be transmitted. Therefore it will be 
difficult to discuss the mode of transmission in term of 
mathematical equations. So for the sake of simplicity, we 
have assumed that 5k is symmetric. Now we define the 
local function / as follows 
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Note that the spatial heterogeneity dispersion changes 
while the local density and individual interaction neigh- 
borhoods arc changed. The discrete random variable 
counts the host in 5k, where Q <nk < 5. Then the mean 
number of hosts per neighborhood is given by 



E[nk{t)] 



N{t)5 
J 



H{t)5- 



(3) 



The rife hosts on 5k include s susceptibles, e exposed 
and i infectives; nk = s -|- e -I- i. Here a is the proba- 
bility, per discrete time interval, that an exposed on 5k 
transmits the disease to the susceptible at site k. /3 is 
the probability, per discrete time interval, that an in- 
fective on 5k transmits the disease to susceptible at site 
k- /i is the probability, per time interval, that an ex- 
posed individuals convert to a infected one. Since the 
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infectives act independently, so the susceptible acquires 
the pathogen in a single time interval with conditional 
probability 1 — (1 — — /?)*. The unconditional prob- 
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ability of infecting the focal susceptible in a discrete time 
interval is e^: 
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rife = n] • Pr [n^ = n] . 
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THE APPROXIMATION OF THE GLOBAL 
DENSITY 



Mean- field analysis have different forms, here we adopt 
a 'local-dispersal' approximation [13]. Let us assume 
S{t), E(t), I(t) are the total susceptible, exposed and 
infected individuals at time t respectively. Then we have 
N{t) = S{t) + E{t) + I{t). The global density of the 
exposed and the infected populations at time t are given 

by 



x{t) 
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(4) 



Then, the global frequency of the exposed and the in- 
fected populations at time t are given by 



E{t) _ E(t) _ E{t) 1 _ x{t) 
N{t) ~ H{t)J ~ ~T~ H{t) ^ H{t)' 



(5) 
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For given independence of sites, the probability that a 
host occupies a given site on 5^ is equal to the global 
density H{t). The individual interaction act as an in- 
dependent event in Sk, so the (e,i) is an independent 
two-dimension variable, and is denoted by (^,77), where 
^ + r] < nk- Form the transmit rule /, we know that 
number of exposed and infected individuals on Sk has a 
(three-parameter) binomial distribution with parameters 
S, x{t), y{t). In the following, the A is the probabihty, 
per time interval, that an infective recovers from infection 
and again become an susceptible one. In the neighbor- 
hood Sk , the exposed and infective at time t+l and t are 
independent. Hence, we can get the global mean density 
of exposed and infective at time t + l are given as fol- 
lows, where x{t + 1) represents the exposed individuals 
and y{t + 1) represents the infected individuals. 



Ijt) _ Ijt) _ m 1 _ y{t) 
N{t) H{t)J J H{t) H{t)' 
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and 



<t + 1) = Fi y{t)) = x{t){l (H{t) - x{t) - y{t)) ■ ^ ^ 

i\) {'~^)x{tfy{t)^(l x{t) - y{t))''^'\l - (1 - a)« . (1 - /?)"] 

5 5 

= xm - + (H(t) x{t) y{t)) ■ tr,^)^(y)^yitr (1 - ^{t) - 

y(i)) -EE© {';,^Hy)^yitr (i - - yit)) - a)« (1 - 

= - /^) + [H{t) - xit) - • [1 - (1 - - /3y(t))' 



(7) 
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yit + l) - F2(^xit),yit)'^ =yit)il-X-d) + ^,xm) 



From Eq. (7) and Eq. (8) we can obtain that the global 
density of exposed and infective are given as follows: 



(9) 



vit + 1)= x{t){l - /i) + i^H{t) - x{t) - y{t) 

1- (l-ax{t) - (3y{t)j 
iy{t + l)=y{t){l-\-d)+iix{t). 
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This is a two-dimension nonlinear system. Since D <C 
N(t) for not too long time, so we can omit the change of 
total individual (see [14]) and let H{t) = H = constant, 
x(t + 1) = x{t) = X, and y{t + 1) = y{t) = y. Where 
the symbol D represents the number of death up to time 
t. The equilibrium points of the system (9) satisfy the 
equations 

X = X ■ {1 - fi) + {H - X -y) ■ [l - {1 - ax - f3yY] , 
y^y{l- j) + fix{t). 

(10) 



Now in the above system , we have 7 = d + A. Since 
6 > 0, < x{t) < H,0 < y{t) < H, so solving these 
above equations we get a zeros-fixed point {xo,yo): 

Xq = 0, 
2/0 = 0. 

Now we consider the Jacobian matrix of the system (9) 
at (0,0). From the equation (7) and the equation (8) 
we obtain 



dFi 
dx{t) 
dFi 

dm 

dF2 
dx{t) 

dF2 
dyit) 



= [l-ax- (3yy + a5{H -x- y){l -ax- PyY'^ - fj. 
= (l-ax- Pyf + I35{H -x- y){l -ax- l3yY~^ - 1, 

= 1-7. 



Therefore the Jacobian Matrix A can be expressed at the 
disease free equilibrium (0,0) as follows 



l + Ha5-n H(i5 
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Then the trace and determinant of the matrix A are given 
by 



trA = 2 - ^ - 7 + H5a 



(11) 



and 



detA = {I- n + H6a){\--i)-H5l3n 

= 1 - n + HSa - -f + fij - HS{aj + /3fi){12) 

Now let us denote TZc = -^^^^^^p^, as the persister 
threshold, then the following theorem holds true. 

Theorem 1. If TZc < 1, then the system (9) is locally 
asymptotically stable around the disease free equilibrium 
(0,0) and unsatble if TZc > 1- 

Proof. Let us linearize the system (9) at (0,0), we get 
jx{t + 1) = (1 + Ha6 - fi)x{t) + HI35y{t), 



y{t + \)^^lx{t) + {\-^)y{t). 



(13) 



Let us consider a positive Lyapunov function 1^ = (^ + 
^)x -f- ^y, then the time derivative of this function along 



the solution of the system (13), we get 



Ay 



V{t + 1)-V{t) 



{-,+^Ht + l)-x{t)\+'-^[y{t+l)-y{t)] 

(3 



/3^ 



Pi 



fl 7 



-)^ix + (- 
7 fj. 



)H5{ax + Py) — (3y + —fix 
7 7 



/3^ 



= -ax - Py+{- + -)H5{ax + Py) 
H 7 

= -{ax + py){l - TZc) < 0. 

Again note that AV = at the origin (0,0). Thus, 
according to Lyapunov LaSallc Theorem [15], the origin 
of system (9) is locally asymptotically stable. 

Now if TZc > 1, one of the Jure conditions [16] is 
violated as 

1 - trA + detA = jfi - H6{aj + Pfi) > 0, (14) 

which implies that the matrix A has a pair of complex- 
conjugate eigenvalues lying outside the unit circle, which 
again implies that the equilibrium (0, 0) of the system (9) 
is unstable. □ 

The above theorem gives us the condition for which 
the pathogen will extinct and when it will persist. The 
system (9) has an infection-free equilibrium in which 
the component of infectives is zero and exposed is zero. 
After analyzing the local stability of the system (9) 
around (0, 0) we get a threshold value which determine 
when the disease dies out and when persist in the 
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system. These threshold conditions are characterized 
by the critical threshold TZc, such that (0,0) is locally 
asymptotically stable liTZc < 1, and unstable if TZc > 1. 
The numerical simulation shows this fact for different 
values of TZr- 



(a) He = 0.4495 



(b) He = 1 




(c) TZc = 7.9653 



(d) Tic = 1.9092 




(e) Tie 



1.0011 



FIG. 4: For each figures we assume the following ini- 
tial condition: a;(0) = 0.001, y{0) = 0.0015. (a) Fig- 
ure a: Figure a have drawn with the parameters val- 
ues: a = 0.002, /3 = 0.004, A* = 0.15,5 = 8,H = 
0.8, d = 0.0003 and A = 0.07. (b) For figure h: a = 
13 = n = 0.64, d = 0.003, and A = 0.637. (c) For 
figure c: with the same parameters value that in fig- 
ure (a) except a = 0.08, (3 = 0.05. (d) For figure d: 
Tic = 1.9092, A* = 0.35, Q = 0.5, /3 = 0.07 and A = 0.45, 
the value of other parameters are remaining same as in 
figure (a), (e) For figure e: TZc = 1.0011,^1 = 0.35, a = 
0.05, /3 = 0.0061105, 2:(0) = 0.01 and y{Q) = 0.015 and 
A = 0.45, the value of other parameters are remaining 
same with figure (a). 



Theorem 2. If TZc = ^'^^7+'^^^ > 1, the system (9) 
possess a positive interior equilibrium {x*,y*). 



Proof. The system of (10) can be expressed as 

Ux ^ (H - X - y)[l - (1 - ax ~ (3yy], 
\y = hx. 



(15) 



Where h — ^/^. Again the equation (15) can be written 



{H ~x- hx) [l - (1 - ax - (3hx)^ 



Now let us assume g{x) = (H — x — hx) [l — (1 — ax — 
l3hx)^] - fix. Since 0<q;<1,0</3<1,0<.t<1 and 
< ?/ < 1, so g{0) = and dg{x)/dx is given by 

g'{x) = -{I + h)[l - {1 - ax ~ (3hx)^] + 

{H -X- hx)6{a + f3h){l ~ ax ~ phx)^'^ - fi. 



Therefore 



5(0) 



H5{a + I3h)- fi 
t^ 

7 



HS{a 
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[HS{a-/ + PfL) - /i7] > 0. 



So, the function g{x) is increasing at (0,0). Therefore if 
we take < r — < 1 , then we get 



9{r) 



< 0. 



Thus, there should be positive interior equilibrium of 
the system in the open interval (0, ^^)- Hence the proof. 



□ 



Now we are interested to find out the positive interior 
equilibrium of the system in open interval (0, jj^)- Now 
for small values of a and /3, taking Taylor expansion and 
keeping only the fist two terms, we get 



(l-ax- Pyf = 1 - 5[ax + f3y) + ■ ■ 
So from the systems (10), we obtain 

\ X ~ x{l — fi) + 5{H — X — y){ax + /3j/), 
1 2/ = 2/(1 - 7) + M^;. 



(16) 



(17) 



Therefore when TZc > 1, we can derive the positive inte- 
rior equilibrium of the system (17) as follows 



Now we are in position to show behaviour of the system 
(9) around the positive interior equilibrium. 



7 M 



7+/^ i5(7+^i)(Q7+/3/i) ' 

Jit 7 

7+M (5(7+/^)(Q7+/3/i.) ■ 



(18) 
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Together with the equation (5) and the equation (6), we 
get 



Solving the equation (24) we get 



_J_M_ 



X 7 

y* _ M 

, H 7+^ H5(7+ai)(q7+/3m) 



ff<5(7+Ai)(Q7+/3/i) ' 

2 

7M 



(19) 



From the equation (19) we can conclude that when 
7?.c > 1, then the global frequency of exposed and infec- 
tive (ff , 7f ) increases with neighborhood size 5 and with 
the pathogen transmission probability a and /3. Figure 
6 shows this fact. Again from Figure 6, we can conclude 
that infections decreases with the increasing of the re- 
covery probability A ( sec Figure 7(a)) but the size of 
the exposed individuals will increase with the increasing 
of A ( see Figure 7(b)) . This observation supports the 
observation of EUner [17]. 

Now Since 0<a<l,0</3<l,0<a;<land 
< 2/ < 1, so < {ax + f3y) < 1. Now if the size of the 
interaction neighborhood and the environment increases, 
that is, (5 ^ J, J ^ oo, then we get 



1-(1 



1, 



as a consequence the equation (10) reduces to 
ix^x{l- fi) + {H -x-y), 

1 y = y(i - 7) + A^a;. 



(20) 



(21) 



Now solving the above equation and denoting the positive 
interior equilibrium as (a;^,?/^) we obtain 



Hi 



M7+7+M ' 



which can be written as follows 



x'^ 7 

H fif+f+fi ' 



(22) 



(23) 



Therefore the frequency of exposed and infected hosts 
are approximated to {[1 + fi+ ^]"\ [1 + 7 + ^]~^)- 
Now the by comparison of the solution (19) with the solu- 
tion (23), it is easy to show that they are approximately 
equal if 5 00 . Note that the solution (19) is better 
than the solution (23), which shows the relation between 
exposed, infected and susceptible individuals. 

Half von Neunann neighborhood: Let us consider 
a special case when S = 2, the neighborhood become one- 
dimensional which is half of the standard von Neunann 
neighborhood. Here in our case the neighborhood have 
only horizontal neighbors or isotropic neighbors, and it 
have the qual probability for each step. Hence we define 
the radius r ~ 1/2, the equation of (10) can be written 
as 



X ~ x{l — /i) + {H — X 



y)[l - {1 - ax - (3yy 



x\ 



yi 



X2 



y2 = 



■y{{2+Ha)j+{2+H/3)iJ.-VA) 

2(7+m)(q7+/3m) ' 
tJ.i{2+Ha)j+{2+Ht3)n-VA) 



■y{{2+Ha)■y+i2+H|3)^^+^/A) 
2(7+M)(a7+/3A') ■ 

IJ.{{2+Ha)-f+{2+HP)fj.+ VA) 
2(7-l-M)(a7+/3M) ' 



where 



A = (iJ/?- 2)2^*2 +72(4- 4iJa- 



H^a^ 



2-tfi{H^al3 - 2H{a + /3) + 2(2 + n)). 



IV. NUMERICAL SIMULATIONS 

In this section we have shown the numerical simulation 
of our CA model with the help Matlab programme. 
Here we have considered the case for a common place 
like any public organization, school, company and so on, 
hence we assume J = 10000. We have also simulated for 
a large value of J but have the similar results like in our 
other paper [18]. The sites are updated synchronously. 
Figure 5 shows that the global frequency of exposed and 
infected individuals increase with the increasing of the 
size of the neighborhood provided TZc > 1 . The Figure 6 
shows that infections decreases with the increasing of 
the recovery probability A ( see Figure 7(a)) but the 
size of the exposed individuals will increase with the 
increasing of A ( see Figure 7(b)) . Figure 7 and Figure 8 
are showing the role of spatial heterogeneity with (5 = 8. 
The Figure 8(a) and figure 8(c) arc showing the disorder 
dispersions of the hosts at initial time. 





(a) 



(b) 



(24) 



FIG. 5: The frequency of exposed and infected indi- 
viduals with a logarithmic scale at mean-field equi- 
librium. The parameters values are: H=0.7, A = 
0.6, d = 0.0003, yU = 0.4 for each figure. In the above 
figures, the upper line is drawn with a = j and 
13 — ^; the middle line with a = i and (3 — j; 
the lowermost line with a = i and P — i. 
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(a) 



(b) 



(a) 



(b) 




FIG. 8: In the above simulations we have taken 
the global density H = 0.7 and the mortality rate 
d = 0.1332 (for Hong Kong city). The value of the 
other parameters are a = 1/8, (3 = 1/5, ^ ~ 0.4, 
A = 0.5. (a) The global frequency of the exposed 
and infected with respect to time, (b) The numbers 
of exposed and infected with respect to time. 



(c) 



(d) 



FIG. 7: In the above figures the space dispersion are 
shown in details with parameters values a — ^, f3 = 
i,Ai = 0.4, d = 0.1332, A = 0.45 and J = 10000. 
Note that a given site may be empty (white), oc- 
cupied by a susceptible (grey), or occupied by an 
exposed, or infected individuals (black), (a) and (b) 
are for low global density H{0) = 0.4, (b) The same 
after 50th iteration of the previous, (c) and (d) are 
for high global density H{0) = 0.85, (d) The same 
after the 50th iteration of the previous. 















— 8^15 

— 8^25 














8^35 
— 8^35 




\ 








































I" 














i" 

































































— 5=15 

— 5=25 











































(a) 



(b) 



FIG. 6: The effects of recovery probability. The 
above figures are drawn for different values of 5 with 
H = 0.7, 11 = 0.8, a = 1/6, d = 0.0003 and /3 1/4. 



In the Figure 7 the white, shaded and black portion 
represent the empty site of susceptible, exposed and 
infected individuals respectively. In Figure 7(a) average 
density of the host is 0.4 and in Figure 7(c) it is 0.85 
at initial time t=0. We have got the Figure 7(b) and 
Figure 7(d) at t=50 with the same parameters values 
as in Figure 7(a) and Figure 7(c) respectively. The 
Figure 8 is the most important figure because it is an 
example of simulation for the SARS disease spreading at 
China in 2003, where the global density is equals to 0.7 
(i.e H = 0.7). The number of the recovery of SARS is 
so small that we can ignore it, so the state of the SARS 
is processing as the Figure 4 shows. The detail rule to 
the CA model is similar to [18]. The Figure 9 is collect 
from the papers of Steven et. al [19] where they have 
studied about the SARS disease in Hong Kong city. Our 
Figure 8 is similar to their figure Figure 9 but not same. 
Similar works have performed by Lipsitch et. al [20]. 
Therefore our results are supported by their works. 
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(a) (b) 



FIG. 9: (a) Distributions for the waiting times of 
the compartments of the stochastic model shown in 
see [19]. Distributions shown are estimates for the 
start of the epidemic. The onset-to-hospitahzation 
distribution changes during the epidemic as a result 
of more rapid hospital admission, (b) Weekly inci- 
dence (by time of hospital admission) with the color 
coding used in (a). 



standard Moore neighborhood. The readers can com- 
pare the another methods of the approximation to spatio- 
temporal models for epidemic with local spread, such as 
PA, SA, HPA in [21]. 

We have made some approximation which are very log- 
ical for the system. The analytical findings and as well as 
the numerical simulations suggest us that the frequency 
of exposed and infected host increase with the increasing 
of heterogeneity (density) and the size of neighborhood 
(See the figure 6, figure 7 and figure 8) . 

It has been observed that to spread the infection, the 
epidemiological threshold TZc > I. This Mean-Field per- 
sister threshold (i?c) property is to some extent related to 
the "basic reproductive ratio i?o" of the epidemic theory, 
see [22]. 

Generally, the effect of dispersion factor depends on 
the size of neighborhood and global density of host. The 
Figure 7 have been performed with different values of 
the global density of the host as the heterogeneity is an 
important factor for epidemic models. 



V. CONCLUSION AND DISCUSSION 

In this article we have proposed and analyzed an epi- 
demic model by using cellular automata theory. We have 
considered here the dispersion of the individuals with in- 
fectious force in latent period. We have also used the 
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